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SIMPLIFIED  MODEL  OF  THERMO-ACOUSTIC  DEVICES 


l 


1  Introduction 


The  analytical  theory  of  thermoacoustic  devices  is  well  developed  for  the  linear  case.  Theoretical  work 
suitable  to  study  higher  amplitudes,  however,  is  in  its  infancy  (Gaitan  and  Atchley  1993). 

To  gain  a  better  semi-quantative  understanding  of  these  processes  as  well  as  of  other  aspects  of  the  design 
of  thermoacoustics  devices  a  simplified  model  suitable  appears  to  be  useful.  The  primary  thrust  of  our  work 
so  far  has  been  the  development  and  study  of  such  a  model.  A  second  aspect  of  our  work  has  been  the 
development  of  the  exact  linear  thoery  along  alternative  lines,  again  with  an  eye  toward  recasting  the  theory 
in  a  different  form  more  suitable  for  extensions  to  the  nonlinear  case. 


Simplified  model  of  thermoacoustic  devices 


The  key  aspect  of  this  model  is  to  recast  the  governing  equations  in  a  form  integrated  over  the  cross-section 
of  the  thermoacoustic  tube  thus  reducing  the  model  to  one  dimension  in  space  (along  the  tube  axis)  and 
time.  Effects  talcing  place  along  the  orthogonal  directions  (friction,  heat  transfer,  etc.)  are  to  be  accounted 
for  by  the  introduction  of  suitable  coefficients  in  the  manner  explained  below. 

Consider  a  linear  thermoacoustic  devices  consisting  of  a  duct  of  variable  area  A(x).  Upon  integrating  the 
equation  of  continuity  over  the  volume  delimited  by  two  neighboring  cross  sections  A(x)  and  A(x  +  dx)  we 
find  the  well-known  form 

dt  a  dx  ’ 

where  p  is  the  gas  density,  u  the  velocity  in  the  x-direction,  and  the  angle  brackets  indicate  the  average  over 
the  cross  section, 


^  A(x)Jm 


dA(.  - .) . 


The  condition  of  that  the  normal  component  of  the  exact  velocity  field  vanishes  on  the  lateral  walls  of  the 
duct  has  been  enforced  in  obtaining  (1). 

Similarly,  the  momentum  equation  in  the  x-direction  becomes 


d  <  pn>  1  dA<  pu2  >  d  <  p>  1  ,  _.  dA 

—w~+a  si  -+~s!r+Ai<’l>-*>-E 


Here  p  is  the  gas  pressure,  r  the  viscous  stress  tensor,  and  the  overline  denotes  the  average  over  the  union 
C  of  the  lines  along  which  the  cross-section  A  is  cut  by  solid  boundaries.  In  the  unobstructed  part  of  the 
duct  this  is  just  the  perimeter  of  A  but,  when  A  is  in  the  stack  region  for  example,  C  would  include  all  the 
“wetted”  perimeter  of  the  stack  cross  section.  Thus,  for  example, 


where  P  is  the  length  of  C  .  The  viscous  component  rxx  has  been  neglected  in  deriving  (3). 

The  last  equation  to  consider  is  that  for  the  gas  enthalpy  that,  with  the  assumption  of  constant  specific 
heat  Cp,  and  the  neglect  of  viscous  heating,  may  be  written  in  conservation  form  as 

UpCpT)  -%  +  pCpTV-u  =  -V-q,  (5) 


where  d/dt  =  d/dt  +  u  •  V  is  the  convective  derivative  and  q  is  the  heat  flux  vector.  Before  averaging  this 
equation  we  note  that,  for  a  perfect  gas, 

pcpt  =  Try?’  («> 
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(7) 


where  7  is  the  adiabatic  index,  so  that  (5)  becomes 

dp 

+  7 pV  u  =  -(7  -  1)V  ■  q . 

Upon  integration  over  the  volume  bounded  by  two  neighboring  cross  sections  we  then  have 

-<p>+-(u.Vp  +  7pV  ■  u)  =  -(7  -  1)  — q  ■  n . 


(8) 


To  make  progress  we  need  to  render  the  previous  expressions  more  explicit.  We  start  by  noting  that  the 
pressure  may  be  taken  to  be  uniform  over  the  cross  section.  This  is  true  in  view  of  the  well-known  properties 
of  viscous  boundary  layers  and  also  because  in  usual  thermoacoustic  devices  the  cross  section  of  the  system 
is  small  compared  with  the  acoustic  wavelength.  With  this  approximation  and  writing  p  in  place  of  <  p  > 
we  have 


(u  ■  Vp  +  7pV  •  u)  ~<  u  >  ^  +  ^pV(A  <  u  >) , 

OX  A 


(9) 


and  also  <  p  >  — p  ~  0.  We  further  disregard  correlation  terms  and  write  the  average  of  products  as 
products  of  averages, thus  dropping  the  angle  bracket  notation.  It  should  be  noted  that  the  error  thus 
incurred  is  greater  than  that  ensuing  from  the  assuming  p  to  be  uniform.  Although  some  correction  could 
presumably  be  introduced  to  mitigate  the  error  we  do  not  do  so  as  our  objective  here  is  to  derive  a  simplified 
mathematical  model.  For  the  same  reason  we  model  the  heat  exchange  between  the  gas  and  the  surrounding 
structure  (including  the  stack)  in  terms  of  a  constant  heat  transfer  coefficient  per  unit  length  H  defined  by 


j<TZ=H(x)[Tw(x)-T\, 


(10) 


where  Tw(x)  is  the  temperature  of  the  solid  structure  (e.g.  the  stack)  at  x.  Ways  to  estimate  H  will  be 
considered  later.  For  the  time  being  note  that  H  can  be  a  function  of  x  to  model  different  intensities  of 
heat  transfer,  e.g.  inside  and  outside  the  stack.  In  a  similar  spirit  we  introduce  fluid-dynamic  resistance 
coefficientper  unit  length  by 

=  -D{x)u.  (11) 


With  the  preceding  simplications  and  approximations  the  system  of  equations  reduces  to 

djP  ,  1  dA'pu  _ 
dt  A  dx 


l+*fc+T75fc4-(T-,)i'tr*w-Tl- 


(12) 

(13) 

(14) 


To  close  this  system,  we  further  postulate  that  the  averaged  quantities  satisfy  the  same  equations  of  state 
as  the  corresponding  exact  ones,  namely 

p  =  RpT,  (15) 

where  R  is  the  universal  gas  constant  divided  by  the  gas  mass. 


3  Energy  relation  and  boundary  conditions 

The  total  energy  of  the  system  £  is  the  sum  of  its  internal  and  kinetic  energies  and,  for  a  perfect  gas,  is 
given  therefore  by 


£  =  f  dx  A(x)  p 
Jx(t) 


CvT+i  u2 


(16) 
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Here  it  has  been  assumed  that  the  left  boundary  oscillates  around  z  =  0  with  and  instantaneous  position  z 
=  X(t).  This  would  be  the  case  if  the  device  is  operated  as  a  refrigerator.  For  prime  mover  operation,  on 
the  other  hand,  X  =  0. 

Upon  taking  the  time  derivative  and  using  the  previous  equations  one  readily  finds 


d£ 

dt 


-  T)  -  Du2}  +  pAX 


x=X(l)  ' 


(17) 


For  a  rigid  termination  at  z  =  L  the  velocity  must  vanish.  As  a  consequence  of  (13),  this  implies  dp'  jdx 
=  0  as  well.  If  Eqs.  (12)  and  (14)  are  then  written  at  z  =  L,  one  finds 

dp  du' 
dt  +  P  dx 


dp 

dt 


+  7  P 


du 

dx 


(7  -  l)/T[T„(z)  -  T] 


z=L 


Upon  eliminating  du/dx  one  then  finds,  at  z  =  L, 


7  p  dT 
7-1  T~dt 


dp 

dt 


+  H(TW  -  T) . 


(19) 


(20) 


This  relation  shows  that  a  knowledge  of  p  at  the  boundary  completely  specifies  T.  No  additional  boundary 
conditions  are  therefore  necessary. 

The  same  argument  applies  at  z  =  0  in  the  prime  mover  case,  which  is  the  only  one  that  we  consider 
from  now  on. 


4  Linear  theory  of  simplified  model 


It  is  useful  to  consider  the  linear  version  of  the  previous  model  in  the  case  of  a  duct  of  constant  area  with 
negligible  dissipation.  In  this  case  Eqs.  (12)  become 


dp'  dppu'  _ 
6t  +  6z  ~ 


(21) 


Here  we  have  written 


du'  dp' 
P°  dt  +  dx 


=  0. 


dp'  du'  .  .  TTrnt 

9r+TW«7  =  -< 


P  —  Po  +  p' .  P  —  Po  +  p' ,  T  —  Tw  +  T1 . 


(22) 

(23) 

(24) 


The  subscripts  0  and  the  prime  indicate  unperturbed  and  perturbed  quantities  respectively.  The  unperturbed 
quantities  satsify  the  equation  of  state  with  po  independent  of  z  so  that  po  =  Po(x)  must  compensate  for  the 
z-dependence  of  Tw . 

It  is  easy  to  eliminate  u'  and  T'  from  these  equations  to  find 


1  \±  (C2  V\  _  1  ,  .  _  n  JL  b 5V  _ 

dt  dx  \  A  dx  J  dt2  +  (7  •  }Rpo  [ 1  dx2 


,d2P'  ay 

dt2 


=  o, 


(25) 


where  c2A  and  c2  are  the  adiabatic  and  isothermal  sound  velocities  given  by 


c2a  =  -yRTw  ,  cj  =  RTW 


(26) 
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Some  insight  into  the  behavior  of  the  solutions  o'  the  linear  problem  can  be  obtained  by  assuming 
oscillations  with  a  (complex)  frequency  ui,  p'  =  P(x)exp  iuit.  With  this  assumption  the  previous  equation 
becomes 


tUl 


'  d 

( ,  QP\ 

1  0  i-» 

,  H 

dx 

\A  $*) 

1  +  W2P 

f'W+uP\ 

=  o. 


(27) 


Upon  multiplication  by  P  (where  the  overline  denotes  the  complex  conjugate),  integration  from  0  to  L,  and 
use  of  the  fact  that  dP/dx  vanishes  at  x  =  0,  L,  we  find  sin  expression  of  the  form 


where 


=  jL\P\ 

Jo 


iu)(Au2  —C)+  Bu ?  —  D  —  0 , 
7-1  fL 


(28) 


'dx 


B  = 


-/*■ 


dP_ 

dx 


dx, 


D  =  - 


j  HTw\P\2dx , 


Po  Jo 
(7-1  )R  fL 


Po 


(29) 

(30) 


It  may  be  noted  that  D  is  in  general  a  complex  quantity  with  real  and  imaginary  parts  given  by 

dx,  (31) 


Jo 


2po 


dP 

2" 

dx 

. 

-=rdP  dP\  , 


(32) 


Equation  (28)  in  general  has  three  complex  roots,  so  that  p'  is  the  superposition  of  three  different  modes. 
For  any  of  these  modes  to  have  the  nature  of  steady  oscillations  it  is  necessary  that  ui  be  purely  real.  In  this 
case,  upon  separating  real  and  imaginary  parts,  we  have  from  (28) 


< j(Acj 2  —  C)  —  Di  =  0 ,  Bui2  —  Dt  =  0  . 


(33) 


It  is  clear  from  the  second  of  these  equations  that  only  two  such  roots  can  exist,  one  the  opposite  of  the 
other.  For  the  first  relation  to  be  satisfied,  it  is  then  necessary  that  £>,•  =  0.  Furthermore,  since  B  is  patently 
positive  definite,  for  the  hypothesized  real  roots  to  exist,  it  is  also  necessary  that  Dr  >  0.  It  is  seen  from  (31) 
that  Dr  is  the  sum  of  a  positive  definite  term  and  a  term  that  can  be  positive  or  negative.  The  requirement 
Dr  >  0  therefore  implies  that  the  latter  term 


J  = 


(34) 


not  be  too  negative. 

When  the  two  equations  (33)  are  compatible,  it  is  readily  shown  that  it  follows  from  (28)  that  the  third 
root  is  purely  imaginary  and  is  given  by 

B 


(35) 


Since,  as  is  clear  from  their  definition,  both  A  and  B  are  positive  definite,  the  mode  corresponding  to  this 
root  is  always  exponentially  damped. 


5  Weakly  nonlinear  theory 

We  now  present  some  preliminary  results  of  a  weakly  nonlinear  analysis  of  the  simplified  model  for  the 
constant-area  case. 
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In  the  first  place,  it  is  useful  to  eliminate  the  density  by  using  the  equation  of  state  to  find 

dp  dT  du  dp  dT 

Ta-”-F  +  pTaI  +  ,‘:r3J-“,,«7  =  0 

du  du  r^srt  dp 

plt+pute+RTfc  =  0 

dp  dp  du 

Furthermore,  by  combining  the  second  and  third  equations  of  this  set  according  to 


d  (9 

RTW  —(Energy)  —  ^-(Momentum) 


we  find  the  wave  equation 


d2u  d2u 
TRT.PgJT-p^T 


dp  ( du  du\  dudu  d2u 

=  m{m+ud^)+pmd^  +  pudIdi  +  R 


?L?1  +  r(t-t  \  d'p 

dtdx  +  R{T  Tw)  dxdt 


du  dp 
dx  dx 


cPp 

‘dx2 


dp  du 
dx  dx 


-RTW  — ^  -  RTvU-rr^  -  yRTw^^  +  (7  -  1  )HRTwi~(Tw  -  T) 


d_ 

dx 


We  now  proceed  in  the  usual  way  introducing  the  slow  times 

r  =  ct ,  .9  =  c2t , 


(36) 


where  f  is  a  measure  of  the  amplitude  of  the  wave,  and  expanding  all  the  fields  in  powers  of  *,  e.g  u  = 
eui  +  c2«2  +  f3«3  +  • . ..  Furthermore,  we  also  let 


(7  -  1)  H  =  h0  +  e/»i  +  e2/»2  +  • . 
Upon  separating  the  various  orders  in  e  we  have: 

Order  e 


(37) 


dpi  dTi  du  1  dTw 

Tw17~p017  +  PoTw  17  - Po17Ul 

_  5ui  .  ffr  9P 1 

po-ft+RTw  — 

dp\  dux 

■9T  +  7P»ar 

a2 ui  a2ui 


=  0 
=  0 
=  0 

—  —RTwho 


dTi 

dx 


Order  c2 

dpi  dTi  du2  dTw 

Twir~P0ir+P()Twi7~P0  17^ 


-T  dPl  I  n  ^ 
Tw  dr  +  Po  dT 


6Pl  m  ,  _  _  .  du\  dpi  dTi  dTw\ 

”  Tllt  +  Pl17  ~ (poTl  +  PlTu>)  17  ~  T”!7Ul  +[P017  +  Pl17)Ul 


3ti2  dp2 

Po  Qt  +  RTW  — 


6 


dn\ ]  dui  dui  t 

>-3T  “Pi “ST  “Pouter  -  RTi- 


if +  7PoS 7 


=  —  {hoTo  +  hiTi)  — —■ 


Ul57-7PlaT 


rvr  w  “2  ' 

jnl  wPQ^T-  ~~  Po ' 


d2U-y 

’  Tl T  ~ 

,  32«l 

2KsfB 

dpi  dui 


D&ui 

yR  dx 2 

<9«i  <9u 


(Tipo  +  Twp\)  +  pi 


1  3<2 

dT\  dpi 


opi  OVl  OUlUll  1  O  Ul  Oil  opi 

+irar + p°-dT'd7 + PoUididi +  R~df  ~di +  m 


prp  d“i  dpi  _  p^r  „  u~pi  -^prr  nr  h  Q±2  j.  h  Ulx 

RTw  dx  dx  1  qx2  7RT*  qx  qx  ^  h°  dx  +hl  dx 


dpi  dui 


)rder  <3 


„  du3  dp2 

Poir  +  RTw  — 


dp3  .  du3 

ar  +  7War 


dui 

du o' 

(  du2 

dux\ 

.  M 

+  ~d7.~ 

{Pl~df 

+  P2~df) 

=  [-  (40t3  +  h ,  t2  +  h ,r, )  -  ^  +  fc)  ]  _  +  u,^  -  7  (p,  |2. + p, 


jKTvPo-z-r  ~  P° 


d2u3 
dt 2 
32ui 


d2ui  32«i  d2u2  92«il  d2u3  d2u  i 

=  2poWfp°a^'  +  2poa^"fPla^t  +Pl~dt^+p2'W 


dpi_dui_  dpi  dui 
dt  dr  dr  dt 


du !  du  1  d2ui  dpi  dTi  d2pi 

P0-5—  -X-  +  P0«ln~a-  +  R-Zr-*Z  +  RT^’aIal 


dx  dr 


dx  dr 


d2pi  (  d2u2 

drdx\  lRTw  VPl  dx2 


dpi  f  du2  du\\  dp2dui  dui  du2  du2dui  du\dui 

+~dt  \jH+Ul!h) +  ~dr~df+podx  dt  +podx  dr+pi~d7~dr 
,  d2**  .  _  d2«i  ,  _  d2“i  .  „9pi  drr2  ,  ndp2dTi 

+PoUl  fcdi +  P0U2dZdi +  PlUldxdt  +R~dI^t+R~dI~dt 

,  „  (r  daP7  d2pi  \  ( du2  dpi  duidpA  f  d2P2  d2pi\ 

+R  \  1  dtdx  +  T2dtdx)  V  dx  dx  +  dx  dx  )  ^  \ 1  dx 2  +  2  dx2  ) 

-*r?T  ( dpi  du 2  -1.  dp2  dui  \  nr  \h  97,3  -l  a  9X2  j.  h  ^ 

-yRTw  V a7  +  ~dZ!7 )  -  ^  r~d7  +  hl-d7  +  h2~dZ 


dx  dx 


The  analysis  of  this  system  is  rather  complex  and  is  still  under  way.  To  demonstrate  the  procedure,  we 
consider  here  the  simplest  case  of  very  weak  heat  transfer. 

Very  weak  heat  transfer 

To  study  this  case  we  set  /i0  =  0,  h\  =  0.  We  start  by  considering  the  Sturm-Liouville  problem 

~Y  +  ^n<r(x)un(x)  =  0 


(38) 


»(0)  =  v{L)  =  0, 

where  the  index  n  labels  the  different  eigenvalues  and  eigenfunctions  and 


Upon  considering  the  problems  for  n  and  m  and  combining  we  have 

Vm  Wn  +  A„<n>„]  -  Un-fam  +  Am<rt;m]  =  0 

(An  Am)  <7Vn  Vm  +  VmVn  VnVm  —  0 
Integrate  over  the  range  the  last  two  terms  to  find 

/  (Vm<  -  «»<)  dx  =  [«„,<  -  Vnv'jfc  -  [  (v'nv'm  -  v'mv'Jdx 

Jo  Jo 

=  [Vmu'n  -  Vndo 

=  0 

The  boundary  conditions  (38)  were  used  in  the  last  step.  Then  we  have 

(An  —  Am)  /  avnvmdx  =  0 
Jo 

which  proves  orthogonality  of  the  eigenfunctions  corresponding  to  distinct  eigenvalues, 


(«n,i>m)H  f  avnvmdx  =  0  if  A„  ^  A„ 
Jo 


Let  us  now  turn  to  the  wave  equation  (38)  to  0{(),  which  is  in  the  present  case 


We  write  the  solution  in  the  form 


d2ui  d2u\ 
,~dx2~.~W 


ui  =  Ai{t,  9)vi(x)  exp(iwt)  +  c.c. , 


where  c.c.  denotes  the  complex  conjugate,  in(x)  is  the  first  eigenvector  of  the  Sturm-Liouville  problem,  and 
A!  =  u,2  the  first  eigenvalue.  Note  that  tq  can  be  taken  real  and  normalized  to  1  according  to  the  scalar 
product  (39).  From  the  linearized  equations  we  then  have 


?  dui 

Pi  =  -7P0-3— +  c.c. 

UI  ox 

dpi  Po  . 

«7  =  + 


92pi  _  po  ( dui  r;  \ 
dx2  UJRTw\dx  TwUl  )+c-c- 

Ti  =  j  (t-IIT.^  +  O,  +  c.c 


where  T. "w  =  dTw/dx. 
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p 


Upon  substitution  of  the  previous  expressions  for  the  first-order  fields,  the  wave  equation  (38)  for  u2 
becomes 

d2u\ 


-yRTu 


d2u2  <32u2 
'dx2  ~  ~dt2~ 


2p0 


drdt 


+tu>p0  ^2(7  -  ljiqvj  +  ~-vt vi'j  [A2  exp(2iut)  -  , 


The  terms  in  the  second  line  drive  oscillations  at  frequency  2w,  which  is  not  an  eigenfrequency  of  the  operator 
in  the  left-hand  side.  Hence,  in  order  to  avoid  secular  terms,  it  is  necessary  and  sufficient  to  impose  that  the 
term  in  brackets  vanish, 


dr 


=  0 


so  that  Ai(r,9)  =  A\(6).  The  method  of  variation  of  constants  enables  us  now  to  determine  the  solution  for 
ui  which  can  be  shown  to  have  the  form 


u2  =  iAi(9)2  u3(x)  exp(uvt)  +  c.c. 


(42) 


where  v2(:e)  is  a  real  function. 

At  the  next  order  we  are  only  interested  in  the  terms  in  the  right-hand  side  of  the  wave  equation  (38) 
that  may  induce  a  resonant  response,  i.e.  the  terms  proportional  to  exp  ±iuit.  Upon  writing  explicitly  only 
such  terms,  the  third-order  wave  equaition  becomes 


yRTu 


d2u3  d2u3 


dx 2 


dt2 


=  2”»"|r '  RT~h'-  dx 


dT, 


—iwpo 


-p o«i 


,  .  ,  ,  du\  9Ty  .  dun 


7(7+5)  /'fc.y  «?  (J  t:  3 rj\  «»,  4-r:c 

5 —  +  s;Us<t+1)+t~  ar“'~5i 


+  c.c. 


9  dA\  t , 
=  2:upoVi— - /i2 

<JU  W 


fRTwrmv\  +  «r.Oi  -  —  ^twvx 


+u)pqA\A\ 

-poA^AiVi 


|(4-3tK 


9T*/„ 


\v 2  ~  yf-viv?  +2viv2 


'2] 


\7(r  +  5)  ,2  r?  [<S_  T»  3T*\  ,4-77; 


(  27./? 


2 


4-  c.c. 


In  order  to  avoid  secular  terms,  the  right-hand  side  must  be  orthogonal  to  tq  in  the  sense  of  the  scalar 
product  (39).  This  condition  readily  leads  to 


dAi  fl2  ~  1  „  .2  7 

^~2a/2pogr^1  2WQi^1 

and  to  the  complex  conjugate  of  this  equation.  Here 
Qr  =  {v\,7RTmT^v\  +  RTwT"v  1  -  - — 


(43) 


g  t> 

Qi  =  w  (  vi,  (4-  37)v,1v2  -  —-vivo  +  2viv2 


7(7  +  5)  ,2 

-  (  vi, - 2 - Ul  + 


w2  T"  3  T'2 


v|_  4-7  r; 

tw  +  2  r. 


v  \ 

fir  /  2  I 

Ul  / 
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Clearly,  the  second  term  of  (43)  describes  a  slight  shift  of  the  linear  frequency.  The  first  term  instead 
indicates  the  possibility  for  damping  ( Qr  <  0),  or  amplification  ( Qr  >  0)  of  the  linear  steady  oscillations. 

6  Numerical  solution 

The  weakly  nonlinear  analysis  previously  described  is  very  cumbersome  and  we  found  it  useful  to  calculate 
some  numerical  solutions  in  parallel  with  it.  For  this  purpose  we  have  taken  H  =  constant  and  we  have 
assumed  a  hyperbolic-tangent  distribution  of  the  wall  temperature, 

Tw  =Tl+  1  +  tanh  1.83178  — ~  —  ~  {T„  -  TL) .  (44) 

Here  the  coordinate  system  used  is  such  that  the  tube  extends  between  -  \L  and  k L ,  with  L=lm.  Figure 
1  shows  the  pressure  (in  atmospheres)  at  the  left-end  of  the  tube  as  a  function  of  the  dimensionless  time 

(45) 

for  L  =  1  m,  H  =  10  W/m  K,  Tl  =  250  K,  T#  =  1400  K.  After  an  initial  transient,  the  result  very  nearly 
stabilizes  around  a  finite  amplitude.  An  example  of  the  pressure  distribution  in  these  steady  state  conditions 
is  shown  in  Fig.  2. 

If  the  higher  temperature  is  increased  to  T»  =  1900 .K,  for  the  same  values  of  H,  L,  and  Tl,  the  resulting 
pressure-versus-time  is  as  shown  in  Fig.  3,  and  the  pressure  profile  in  the  tube  at  a  particular  instant  as 
shown  in  Fig.  4.  The  oscillations  exhibited  in  this  case  are  clearly  a  numerical  artifact  due  to  the  use  of 
upwind  differencing  in  the  spatial  discretization  of  the  equations.  This  problem  is  well  known  to  occur  in 
the  neighborhood  of  shock  fronts.  A  new  class  of  numerical  methods,  called  total-variation-diminishing,  has 
been  developed  to  avoid  it.  We  are  currently  implementing  one  such  technique  to  be  able  to  extend  the 
calculation  to  higher-amplitude  oscillations. 

Reference 


Gaitan,  D.F.  and  Atchley,  A. A.  1993  Finite  amplitude  standing  waves  in  harmonic  and  anharmonic 
tubes.  J.  Acousi.  Soc.  Am.  93  2489-2495. 
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Fig.  1.  Pressure  at  the  left  wall  of  .the  tube  (in  atmospheres)  versus  the  nomlimensional 
time  defined  by  (45)  for  the  temperature  distribution  given  by  Eq.  (44)  witli  Tl  =  250  K, 


Fig.  2.  Pressure  distribution  (in  atmospheres)  Tor  the  case  of  the  previous  figure  at  the 
dimensionless  time  t .  =  238.5. 


C.  Herman  and  C.  Bartscher 


EXPERIMENTAL  EFFORT 


1.  INTRODUCTION 


The  main  objective  of  the  experimental  work  is  to  uncover  the  factors  limiting  the 
performance  of  existing  designs  of  thermoacoustic  devices  by  gaining  better  insight  into 
the  flow  structure  and  the  heat  transfer  mechanisms.  This  will  be  achieved  by  using 
different  visualization  techniques,  such  as  smoke  injection,  holographic  interferometry  and 
thermochromatic  liquid  crystals.  Visualization  techniques  will  provide  the  necessary 
information  on  the  local  flow  and  temperature  fields  in  the  regions  of  interest  in  addition 
to  the  conventional  measurements  of  temperatures,  pressures  and  other  operating 
parameters  in  different  parts  of  the  thermoacoustic  device.  The  measurements  to  be 
performed  during  the  present  study  differ  from  those  described  in  literature,  and  the 
requirements  imposed  by  the  different  visualization  measurements  require  specific  design 
solutions  for  the  experimental  setup.  The  specific  design  requirements  and  procedures  are 
described  in  this  report 

2.  DESIGN  OF  THE  EXPERIMENTAL  SETUP 

The  design  tasks  for  the  setups  to  be  used  in  the  experimental  investigation  of 
thermoacoustic  phenomena  are  twofold: 

1,  In  the  first  step  of  the  experimental  investigations,  an  enlarged  model  of  a 
thermoacoustic  device,  suitable  for  visualization  of  the  oscillatory  temperature  fields 
by  holographic  interferometry,  was  built  The  model  was  designed  so  that  with  minor 
modifications,  it  can  function  both  as  a  heat  pump  and  a  prime  mover.  The  design  also 
allows  the  visualization  of  the  flow  fields  by  smoke  injection  and  the  visualization  of 
the  temperature  distribution  on  the  stack  plates  using  thermochromatic  liquid  crystals. 

2.  At  a  more  advanced  stage  of  the  research  (second  year),  after  satisfactory  results  have 
been  obtained  in  the  experiments  with  the  enlarged  model,  a  real  thermoacoustic 
device  suitable  for  accurate  measurements  will  be  built.  The  design  should  be  modular 
in  order  to  allow  the  modification  and  exchange  of  its  key  parts.  In  this  way,  it  will  be 
possible  to  investigate  modifications  proposed  during  the  course  of  the  research  on  the 
actual  device.  Again,  an  attempt  will  be  made  to  build  a  device  which  can  function 
both  as  a  prime  mover  and  a  heat  pump.  This  device  will  be  suitable  for  the 
investigation  of  the  phenomenon  of  thermoacoustic  streaming. 

In  addition  to  the  two  experimental  setups,  a  thermoacoustic  couple,  suitable  for 
temperature  measurements  both  in  the  thermoacoustic  device  and  in  the  model,  will  be 
designed.  Measurements  and  visualization  experiments  with  the  thermoacoustic  couple 
only  will  be  performed  to  analyze  heat  transfer  and  fluid  flow  phenomena  on  a  short  stack. 
The  thermoacoustic  couple  will  be  designed  to  allow  visualization  experiments  with  smoke 
injection  as  well  as  interferometric  measurements. 
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3.  ENLARGED  MODEL  OF  A  THERMOACOUSTIC  DEVICE 


3.1  Design  requirements 

Both  standard  requirements  related  to  the  operation  of  thermoacoustic  devices  and 
specific  ones,  which  allow  visualization  measurements,  have  to  be  met  in  the  design  of  the 
experimental  setup.  The  most  important  criteria  for  the  design  of  the  enlarged 
thermoacoustic  device  model  are  as  follows. 

1.  A  modular  structure  allows  the  experimenter  to  disassemble  the  device  completely 
and  easily  in  order  to  exchange  its  key  parts.  In  this  way,  the  flexibility  to  investigate  a 
variety  of  heat  transfer  and  fluid  flow  phenomena  under  different  operating  conditions 
is  achieved.  At  this  stage,  we  must  keep  in  mind  that  the  need  to  investigate  new 
effects  may  arise  in  different  phases  of  the  research. 

2.  In  order  to  allow  measurements  by  holographic  interferometry,  the  flow  and 
temperature  fields  in  the  investigated  region  should  be  two-dimensional.  This  is 
achieved  by  designing  the  experimental  setup  in  such  a  way,  that  the  cross  section  of 
the  resonance  tube  is  rectangular.  However,  all  thermoacoustic  devices  previously 
described  in  literature  have  a  circular  cross  section.  The  temperature  in  the  tube  is 
constant  along  the  paths  of  the  laser  light  used  for  transillumination  of  the 
measurement  volume.  The  requirements  for  holographic  interferometry  define  the 
spanwise  dimension  of  the  tube.  For  the  temperature  differences  expected  in  the 
experiments,  a  spanwise  dimension  of  15  cm  is  required.  The  height  of  the  tube  is 
determined  by  the  maximum  efficient  diameter  of  the  expanded  laser  beam,  which  is 
about  5  cm.  Using  the  above  criteria,  the  dimensions  of  the  cross  section  of  the 
resonance  tube  are  selected  to  be 

•  height  of  tube:  h  0.05  m 

•  width  of  tube:  w  0.15  m 

3.  The  plates  of  the  stack  are  parallel  in  order  to  achieve  approximately  two- 
dimensional  temperature  fields  in  the  region  of  interest.  The  minimum  spacing  between 
the  stack  plates  is  determined  by  the  resolution  of  holographic  interferometry.  In  our 
investigations,  this  spacing  is  determined  to  be  5  ram. 

4.  The  viewing  windows  allow  the  visualization  of  a  region  of  8  cm  in  diameter,  which 
corresponds  to  the  diameter  of  the  expanded  laser  beam.  The  region  of  interest  for 
visualization  of  flow  and  temperature  fields  includes  the  stack  and  the  channel  region 
in  the  vicinity  of  the  stack,  and  is  approximately  12  cm  long.  As  the  dimension  of  this 
region  is  larger  than  the  diameter  of  the  expanded  laser  beam,  one  half  of  this  region  is 
investigated  at  a  time.  Two  sets  of  interferometric  images  and  the  corresponding 
viewing  windows  are  necessary  in  order  to  record  temperature  fields  along  the  whole 
stack  and  in  its  vicinity.  The  direction  of  transillumination  is  indicated  in  the  schematic 
of  the  experimental  setup  in  Figure  1. 

5.  The  tube  walls  in  the  region  of  the  stack  and  the  stack  plates  are  transparent  to 
allow  flow  visualization  by  smoke  injection  as  well  as  the  measurement  of 
temperatures  of  the  stack  plates  by  coating  them  with  thermochromatic  liquid  crystals. 
In  both  of  these  measurements,  the  direction  of  transillumination  is  perpendicular  to 
that  during  interferometric  measurements,  as  indicated  in  Figure  1. 
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6.  The  channel  will  be  built  with  a  closed  end.  and  with  constant  cross  section  of  the 
resonance  tube.  The  selected  design  allows  the  variation  of  the  tube  length,  if 
necessary. 

7.  The  model  operates  with  air  at  atmospheric  pressure  as  the  working  fluid.  For  the 
design  calculations,  the  following  thermophysical  properties  of  air  are  assumed 

•  mean  temperature:  Tm  293  K 

•  ratio  of  specific  heat:  y  1.40 

•  specific  gas  constant:  R  286.6  J/Kg  K 

•  velocity  of  sound:  cUT  342.9  m/s 

8.  Apart  from  the  visualization  measurements,  pressure  measurements  (to  obtain 
information  on  losses  in  the  acoustic  field  due  to  the  stack  and  the  heat  exchangers), 
temperature  measurements  -  with  thermocouples  -  along  the  surface  of  the  tube  (to 
account  for  heat  losses),  temperature  measurements  in  the  heat  exchangers  and  near 
the  acoustic  driver,  and  other  measurements  are  necessary. 


Figure  1.  Schematic  of  the  thermoacoustic  device  model  with  viewing  directions  for 
visualization  measurements 

3.2  Acoustic  driver 

A  commercial  loudspeaker,  an  8"  woofer,  is  used  as  acoustic  driver  (Peerless  -  type 
832556).  The  most  important  technical  data  for  the  selected  loudspeaker  are  the 
following: 

•  frequency  range:  /  20  -  5000  Hz 

•  resonant  frequency:  F0  30.9/32.3  Hz 

•  effective  cone  diameter  d  0.153  m 

•  max.  linear  SPL:  108  dB/85  W 
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33  Calculation  of  the  length  of  the  resonance  tube  and  the  frequency  range 
With  the  assumptions  discussed  in  the  previous  sections,  it  is  possible  to  perform  the 
necessary  calculations  to  determine  the  limiting  dimensions  of  the  resonance  tube  and  to 
verify  that  the  selected  loudspeaker  can  deliver  the  required  acoustic  power  input. 


Minimum  and  maximum  tube  length 

a)  Plane  acoustic  wave  condition.  For  better  efficiency,  plane  acoustic  wave  fronts 
are  preferred  during  the  operation  of  the  thermoacoustic  device.  In  order  to  achieve  this, 
the  diameter  of  the  tube  must  be  sufficiently  small  compared  to  the  dimensions  of  the 
loudspeaker,  and  the  ratio  of  the  wavelength  of  the  sound  wave  to  the  diameter  of  the  tube 
must  be  greater  than  six  [1].  As  the  cross  section  of  the  resonance  tube  is  rectangular,  the 
larger  dimension  (width  of  the  tube,  w=0.15  m)  is  selected  as  representative  value  for  the 
diameter.  For  this  situation,  the  minimum  wavelength  is  determined  to  be 


—  >  6 
w 


=  0 .  90  m. 


b)  Maximum  frequency  for  a  flat  reproduction.  The  loudspeaker  has  perfect 
compensation,  and  achieves  a  "flat"  response  up  to  a  specific  maximum  frequency.  The 
following  equation  shows,  that  an  increase  of  the  loudspeaker  diameter  d  results  in  a 
decrease  of  the  available  maximum  frequency  for  the  flat  reproduction  [2]. 


f tnx 


4n±d 


342.9? 
4rt40.153m 


=  356.69  Hz 


Thus,  the  upper  frequency  limit  in  order  to  achieve  a  flat  reproduction  is  356  Hz  for  the 
selected  loudspeaker  with  a  diameter  d=0. 153  m.  This  calculation  holds  for  a  temperature 
of  293  K  only.  As  the  velocity  of  sound  is  dependent  on  the  temperature,  the  value  for  the 
frequency  will  also  change. 

c)  Minimum  tube  length.  From  the  maximum  frequency,  the  lower  limit  for  the  tube 
length  l  is  determined  at  room  temperature  from  the  following  equations: 


342.9? 


=  0.96  m  and  /. 


tube.min 


=  0 .  48  m 


/«  356.69 Hz 

d)  Maximum  tube  length.  As  the  resonant  frequency  for  speaker-type  832556  is 
32.9  Hz,  the  minimum  frequency/^  =  35  Hz  is  used  in  the  calculation. 


c.  343 11* 

=  -^g-=  '  =9.80  m 

f  ■  35 1 

J  min  s 


and  /, 


tube, max 


=  4 .  90  m 


Length  of  the  resonance  tube  and  the  frequency  range 

As  the  maximum  possible  tube  length  of  4.90  m  is  not  suitable  for  interferometric 
measurements,  the  following  dimensions  of  the  tube  length  are  selected  for  our 
measurements 

0.5334  m  <  /  <  0.9144  m. 

Hence,  the  range  of  operating  frequencies  is 

187  Hz  <  /  <  321Hz. 


16 


3.4.  Optimum  stack  position  and  expected  temperature  differences 
Theory 

Atchley  et  al.  (1990)  have  conducted  a  thorough,  quantitative  investigation  of  the  basic 
theory  underlying  thermoacoustic  heat  transport  for  the  simplest  class  of  thermoacoustic 
engine,  a  stack  of  short  plates  referred  to  as  a  thermoacoustic  couple  (TAC).  They 
measured  the  temperature  difference  developed  across  the  short  stack  of  plates  in  acoustic 
standing  waves  as  a  function  of  the  position  of  the  plates  in  the  acoustic  field,  drive  ratio, 
plate  configuration,  thermal  properties  of  the  plates,  and  thermophysical  properties  of  the 
gas.  The  temperature  differences  between  the  two  ends  of  the  stack  plates  are  due  to  the 
net  effect  of  the  shuttling  action  of  the  gas  parcels  by  removing  heat  from  one  side  and 
accumulating  the  heat  at  the  other.  In  other  words,  the  end  of  the  plate  closer  to  the 
pressure  node  gets  cooler  and  the  side  of  the  stack  closer  to  the  pressure  antinode  gets 
warmer  and  there  is  no  net  heat  loss  or  gain  by  interior  portions  of  the  plates. 

The  equation  for  the  theoretical  temperature  gradient 

The  following  equation  for  the  steady-state  temperature  difference  developed  across  a 
TAC  in  an  acoustic  standing  wave  was  originally  derived  by  Wheatly  et  al.  (1983)  and 
slightly  altered  by  Atchley  et  al.  (1990)  to  include  the  effect  of  thermal  conduction 
through  the  gas: 


A  T  = 


1 


K. 


*  Ax*sin(47cX  ‘x) 


p  c  K-  4P mc(Kpdp  +  Kgdg)(l  +  o) 

r  mp  ^ 


1+ 


K. 


P02(l-aVa) 


Pmcpn±  4( Kpdp  +  Kgdg  )p  Jm 2n^-(y- 1)(1  -  a2 ) 


[l-cos(47tA.~1x)] 


For  this  equation  it  is  assumed  that  Ax«X/2 n;  8K  ,  5V  «  dg ;  and  AT  «  Tm  and  that  the 

presence  of  the  stack  does  not  significantly  influence  the  acoustic  standing  wave. 
Furthermore,  the  heat  capacity  of  the  plate  is  considered  much  greater  than  that  of  the  gas, 
so  that  the  temperature  of  the  plate  surface  does  not  oscillate  on  an  acoustic  time  scale. 

In  order  to  estimate  a  suitable  stack  length  and  stack  position  as  well  as  to 
calculate  the  temperature  difference  corresponding  to  different  stack  length  and  stack 


positions  following  variables  and  parameters  are  determined: 
a)  Variables 

AT  temperature  difference 

°C 

result 

Ax 

stack  length 

m 

0.01  -0.10  m 

X 

position  of  center  point  of  stack  from  rigid  wall 

m 

0- 0.0254  1 

Po 

dynamic  pressure  amplitude 

Pa 

0.5%  -  3  % 

1 

tube  length 

inch 

21  to  36 

Tm 

mean  temperature 

K 

273  K;  293  K 
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b) 

Constants 

c 

speed  of  sound  in  air 

m/s 

depends  on  T, 

Cp 

isobaric  heat  capacity  per  unit  mass  for  air 

kg/m3 

1006 

dg 

plate  separation 

m 

0.005 

dp 

plate  thickness 

m 

0.001 

X 

wavelength 

m 

0.0508  •  1 

7 

ratio  of  specific  heats  of  air 

- 

1.40 

kg 

thermal  conductivity  of  air 

W/mK 

depends  on  T, 

kp 

thermal  conductivity  of  plate  (acrylglass) 

W/mK 

0.184 

Pm 

mean  air  density 

kg/m3 

depends  on  T, 

a 

Prandtl  number  of  air 

- 

0.72 

Temperature  differences  for  different  stack  lengths 

Figure  2  presents  graphically  the  result  of  the  calculation  of  the  theoretical  temperature 
differences  by  varying  the  stack  lengths  between  0.01m  and  0.1m,  at  a  frequency  of  321 
Hz  and  corresponding  tube  length  of  0.5334  m.  The  operating  temperature  selected  for 
the  calculations  is  293  K,  which  means  that  the  speed  of  sound  in  air  is  c  =  342.9  m/s,  the 
thermal  conductivity  of  air  is  Arg  =  0.0257  W/mK,  and  the  mean  density  of  air  is  pm  =  1.205 
kg/m3.  The  maximum  dynamic  pressure  amplitude  is  considered  to  be  3%  of  the  mean 
atmospheric  pressure. 

The  calculations  indicated  that  by  increasing  the  stack  length  the  theoretical 
temperature  difference  increases  as  well.  In  order  to  satisfy  both,  a  maximum  theoretical 
temperature  difference  and  the  requirements  for  the  interferometric  measurement 


Figure  2:  Theoretical  temperature  difference  by  varying  the  stack  length  and  stack 
position  within  the  resonance  tube 


18 


Determination  of  the  exact  stack  position 

The  optimum  position  of  the  stack  should  be  a  compromise  to  match  the  requirements 
regarding  different  temperatures,  pressure  amplitudes,  tube  lengths  and  frequencies.  In 
order  to  satisfy  these  requirements  several  cases,  applying  the  equation  mentioned  above, 
are  calculated.  The  best  compromise  could  be  achieved  by  a  stack  center  position  of  0. 127 
m  from  the  loudspeaker  membrane,  which  ensures  high  enough  temperature  differences 
between  the  ends  of  the  stack  plates. 


4.  INSTRUMENTATION 

Acoustic  power  delivered  by  the  loudspeaker 

For  calculations  of  the  acoustic  power,  the  knowledge  of  the  oscillation  amplitude  of  the 
loudspeaker  membrane  is  essential.  The  minimum  amplitude  is  obviously  0  mm,  and  the 
maximum  amplitude,  indicated  by  the  loudspeaker  manufacturer,  is  6  mm.  For  the 
estimation  of  the  theoretical  acoustic  power  input  to  the  resonance  tube,  the  appropriate 
working  range  for  the  membrane  amplitude  was  selected  as: 

amplitude:  A^  =  1.0  mm  and  A^  =  4.0  mm, 

and  the  limiting  frequencies  are 

minimum  frequency  =  187  Hz 

maximum  frequency  =  32 1  Hz 

With  the  following  equations 

•  s(t)  =  A  sin(tor) 

•  vft)  =  A  to  cos(cor) 

•  aft)  =  -  A  to2  sin(cor)  and:  a)  =  constant  =  2nf 

with:  A  -  amplitude 

to-  angular  frequency 

both  velocity  and  acceleration  of  the  loudspeaker  membrane  can  be  calculated. 

Velocity  of  the  loudspeaker  membrane 

The  maximum  value  for  the  velocity,  v^,  is  determined  from  the  condition  cos(cor)  =  1. 
Hence,  the  time  t  is: 

t  =  —T  with:  n  =  1,2,3, . 

2 

Thus,  the  maximum  velocity  of  the  loudspeaker  membrane  will  be: 

vm»=ACO=2A7t/, 

which  leads  to  the  following  results  for  the  maximum  velocity  for  the  given  ranges  of 
frequencies  and  amplitudes: 

=>1.175  -?-  <  <  4.70  -f-  for /=  187  Hz  and  0.001  m  <  A  <  0.004  m 

=>2.017  <  8.06  •*-  for /=  321  Hz  and  0.001  m  <  A  <  0.004  m 

Acceleration  of  the  loudspeaker  membrane 

In  order  to  obtain  the  maximum  value  for  the  acceleration  of  the  loudspeaker  membrane 
amw’  the  condition  sin(cor)  =  1  has  to  be  satisfied.  Hence  the  time  t  is 
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t  =  - -  T  with  n=  1,2,3. 

4 

The  maximum  acceleration  will  be 

<2miX  =  /4  (D2  =  4  A  n2P 


=>  1380  -4-  <  <  5522  for /=  187  Hz  and  1.0  mm  <  A  <  4.0  mm 

=>  4068  <  16271  -f-  for /=  321  Hz  and  1.0  mm  <  A  <  4.0  mm 


Acoustic  energy  generated  by  the  loudspeaker 

The  acoustic  energy  that  can  be  theoretically  radiated  into  the  tube  of  the  thermoacoustic 
heat  pump  model  is  calculated  as 

P„  =  ->r  nr2pcR{U2, 


with 


4 


(2  kr)2 
8 


The  values  used  in  the  calculations  are  the  following: 

radius  of  the  loudspeaker  membrane  r  =  0.076  m 
characteristic  impedance  of  air  [5]  pc  =414.0-^- 


speed  of  sound 
range  of  frequencies: 

Hence: 


c  =  343.11 

187  Hz  <  /  <  321  Hz 


range  of  variable  Rx  with  k  =  -f- : 
range  of  membrane  velocities: 


0.03386674  <  R{  <  0.099793 
1.175  -f  <U<  8.06  -f- 


For  these  values  the  power  of  the  sound  emitted  by  the  membrane  is 

0.176  <PK<  24.35  W 


Expected  pressure  amplitudes 

The  next  step  in  the  calculation  is  to  estimate  the  pressure  amplitude  that  can  be  delivered 
from  the  membrane  of  the  loudspeaker.  The  pressure  can  be  determined  as 

p  =  v  pc 

with:  p  -  pressure  difference  relative  to  atmospheric  pressure 

v  -  membrane  velocity  and 
pc-  characteristic  impedance  of  air 

As  the  most  important  values  for  this  calculation  are  the  minimum  and  maximum 
membrane  velocities,  =  1.175-“-  and  =  8.06  -f-.  The  pressure  amplitude  to  be 

expected  is  in  the  following  range: 

486  Pa  <  p  <  3337  Pa 
The  dynamic  pressure  is  then  in  a  range  of 

0.48%  -  3.3% 

relative  to  the  atmospheric  pressure  of  pm  -  101325  Pa.  These  calculations  indicate  that 
the  selected  loudspeaker  is  capable  of  delivering  sufficient  acoustic  energy  to  drive  the 
thermoacoustic  heat  pump  model. 
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Measurement  of  acoustic  energy 

The  pressure  transducer.  In  order  to  calculate  the  actual  acoustic  power  input  during  the 
operation  of  the  heat  pump,  information  on  the  dynamic  pressure  amplitude  is  required. 
Pressure  transducers  use  sensor  chips  or  fragile  membranes  to  measure  pressure. 
Differential  pressure  transducers  measure  the  difference  between  two  pressure  sources. 
The  pressure  difference  is  then  obtained  by  pointing  one  side  of  the  device  membrane  to 
one  pressure  source  and  leaving  the  other  port  open  to  the  atmosphere.  As  the  maximum 
change  of  the  dynamic  pressure  is  about  3%  of  the  atmospheric  pressure,  a  gage  pressure 
transducer  is  the  obvious  choice. 

One  pressure  transducer  that  fulfills  the  requirements  regarding  accuracy,  cost,  low 
noise,  air  flow  monitoring  and  low  power  consumption  is  the  SX01D  pressure  transducer 
from  SenSym.  This  transducer  is  operating  in  a  differential  pressure  range  between  0  and  1 
psi.  This  matches  well  the  operating  pressure  range  of  the  heat  pump. 

The  accelerometer.  In  order  to  determine  the  sound  power  emitted  by  the  acoustic  driver 
during  the  operation  of  the  heat  pump,  the  velocity  of  the  loudspeaker  membrane  must  be 
measured.  As  the  velocity  cannot  be  measured  directly,  an  accelerometer  will  be  used,  and 
the  acceleration  signal  will  be  integrated  into  velocity.  An  important  requirement  is  a  low 
mass  of  the  transducer,  in  order  to  avoid  influencing  the  membrane  mass. 

The  accelerometer  chosen  for  the  present  application  is  the  EGA- 125-2500 
Miniature  Accelerometer  from  Entran.  Its  weight  without  wiring  is  0.5  grams,  it  has  an 
acceleration  range  of  ±  2500g  with  a  sensitivity  of  0.069  mV/g  and  can  be  used  up  to  a 
frequency  of  1800  Hz.  These  specifications  exactly  match  all  the  requirements  for  our 
application.  The  semiconductor  circuitry  of  the  transducer  is  fully  compensated  for 
temperature  changes  within  a  range  from  21°  to  77°C. 

The  power  supply.  In  order  to  drive  both  the  pressure  transducer  and  the  accelerometer, 
the  PS- 15  power  supply  by  Entran  was  selected.  It  has  an  adjustable  sensor  excitation 
from  1  to  15  VDC,  and  can  be  used  for  both  transducers.  The  power  supply  generates  a 
DC  supply  voltage  of  high  stability  and  low  noise.  The  signal  output  port  of  the  power 
supply  will  mainly  be  used  for  the  accelerometer,  whereas  the  signal  output  from  the 
pressure  transducer  will  be  connected  to  an  amplifier  or  a  data  acquisition  board. 

5.  ACCOMPLISHMENTS 

During  the  past  project  period,  the  thermoacoustic  device  model  was  built,  measurement 
tasks  were  identified,  and  sensors  necessary  for  the  measurements  were  selected  and 
purchased.  A  photograph  of  the  thermoacoustic  device  model  and  instrumentation  is 
presented  in  Figure  3. 
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Figure  3.  Thermoacoustic  device  model  and  instrumentation 


6.  FUTURE  PLANS 

Experiments  with  the  thermoacoustic  device  model: 

1.  Basic  power,  pressure  and  temperature  measurements. 

2.  Visualization  of  temperature  fields  by  holographic  interferometry  and  heat  transfer 
measurements. 

3.  Flow  visualization  by  smoke  injection. 

4.  Measurement  of  the  temperature  distribution  in  the  stack  plates  using  thermochromatic 
liquid  crystals. 

5.  Evaluation  of  the  visualization  measurements. 

6.  Design  of  the  heat  exchangers. 

7.  Investigation  of  the  operating  mode  as  prime  mover. 
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OBJECTIVES 

The  objective  of  this  phase  of  the  effort  is  the  development  of  a  computational  methodology 
for  simulation  of  the  flow  in  the  neighborhood  of  a  thermoacoustic  stack,  and  the  implementation 
of  the  numerical  schemes  to  analyze  flowfield  response  to  different  acoustic  forcing  levels. 

This  effort  complements  ongoing  experimental  and  theoretical  efforts  aiming  at  an  improved 
fundamental  understanding  of  thermoacoustic  devices.  The  development  of  these  parallel  efforts  is 
described  separately.  This  section  thus  provides  a  summary  of  numerical  modeling  activities. 


APPROACH 

In  its  simplest  form,  a  thermoacoustic  devices  consists  of  a  straight  resonance  tube  and  a  stack 
of  parallel  plates.  The  arrangement  is  schematically  illustrated  in  Fig.  1.  By  exciting  a  standing 
wave  within  the  resonance  tube,  it  is  now  well-established  that  a  temperature  gradient  develops 
across  the  stack,  thereby  enabling  heat  exchange  between  a  two  systems  or  between  a  system  and  a 
reservoir. 


Figure  1.  Schematic  illustration  of  a  thermoacoustic  device. 


As  mentioned  above,  simulation  of  the  flowfield  in  a  small  neighborhood  of  the  stack  is 
targeted.  Since  the  length  of  the  stack,  L ,  is  typically  much  smaller  than  that  of  the  resonance  tube, 
A,  it  is  highly  desirable  to  adopt  a  simplified  flow  model  which  enables  us  to  isolate  the  stack, 
while  at  the  same  time  retaining  all  of  the  essential  features  of  the  entire  system.  To  this  end,  we 
first  adapt  a  simplified  Iow-Mach-number  compressible  model  which  meets  this  requirement 


Figure  2.  Schematic  illustration  of  a  thermoacoustic  stack  consisting  of  thin 

closely  spaced  flat  plates. 
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The  basic  assumptions  used  in  the  construction  of  the  model  are:  (1)  the  length  of  the  stack  is 
much  shorter  than  that  of  the  resonance  tube  so  that,  in  a  small  neighborhood  of  the  stack,  the  low- 
Mach-number  limit  of  the  compressible  momentum  and  energy  equations  can  be  adopted;  (2)  away 
from  the  stack,  the  flowfield  is  assumed  to  be  well  approximated  by  an  idealized  acoustic  standing 
wave.  These  assumptions  enable  us  to  follow  a  construction  that  is  similar  to  that  used  in  the 
development  of  a  large  number  of  reacting  flow  models,  e.g.  [1-4].  We  start  from  the  governing 
equations  for  compressible  non-reacting  flow  of  an  ideal  gas  (Table  1),  normalized  with  respect  to 
appropriate  characteristic  length,  velocity,  density  and  pressure  scales. 


TABLE  1 

Governing  Equations  for  Compressible 
Fluid  Flow 


Mass 

+ pV  •  u  =  0 

Dt 

(1.1) 

Momentum 

yM*  =_V p+>^1V.T 

1  H  Dt  y  Re 

(1-2) 

Energy 

DT  y-IDp  1  _  (y-l)M2 

H  Dt  y  Dt  Pe  H  Re  Y 

(1.3) 

State 

•a 

II 

(1.4) 

Next,  we  expand  all  relevant  quantities  in  terms  of  powers  of  e  s  yM1,  where  y  is  the  specific 
heat  ratio  and  M  is  the  Mach  number,  i.e.  by  expressing  a  generic  gasdynamic  variable  £(x,f)  in 
series  form  as: 


S(x,r)  =  qo  (X,  t)  +  eqx  (x,0  +  £2£2(x,f) + •  •  • 


(1) 


Substituting  the  appropriate  expansions  into  the  momentum  equations,  and  collecting  the  zeroth 
power  in  £,  immediately  yields; 


Vp„  =0 


(2) 
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i.e.  the  zeroth  component  of  pressure  is  spatially  uniform,  and  therefore  a  function  of  time  only. 
By  an  abuse  of  terminology  and  notation,  we  shall  call  this  pressure  component  die  thermodynamic 
pressure  and  denote  it  by  /*(/).  To  first  order  in  £,  one  has: 


Momentum  (e): 


+  — V*T 
Re 


(3) 


Energy  (£): 


Izi^  =  JLv’r 

y  dt  Pe 


(4) 


Note  that  the  energy  balance  is  independent  of  pl  and  that  viscous  dissipation  has  been  neglected. 
A  more  convenient  form  of  the  energy  equation  is  obtained  by  differentiating  the  (zeroth  order) 
equation  of  state: 


dP 

dt 


DTa  „  Dp 
P,-^  +  T. 


Dt 


Dt 


(5) 


and  combining  with  the  continuity  equation  to  get: 


V 


•  u  =  — 

u<» 


yP  dt+  PPe 


v2r 


(6) 


In  the  geometry  considered,  velocity  fluxes  into  the  system  are  known,  following  the  assumption 
that  matching  with  an  idealized  sound  wave  away  from  a  stack  occurs.  This  enables  us  to  integrate 
the  above  equation  over  the  entire  volume  of  the  computational  domain  to  get: 


~—  =  -P\u*ndA  +  —[VT*ndA 
y  dt  J  PeJ 


(7) 


Next,  we  adopt  a  vorticity-based  formulation  of  the  equations  of  motion.  The  derivation  relies  on 
the  Helmholtz  decomposition  of  the  velocity  field  into  divergence-free  and  irrotational  components, 
and  on  the  vorticity  form  of  the  momentum  equation.  Thus,  we  first  let: 


u„  =  +  u# 


(8) 


where  ^  is  a  potential  function  and  uM  is  the  divergence-free  vortical  velocity  component. 
Substituting  Eq.  (8)  into  Eq.  (6),  one  gets: 
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Second,  taking  the  curl  of  the  momentum  equations,  we  have: 


^■  +  Vx(ax  u)  =  -^x-t— Vx(Vu)  (10) 

3,  p.  Dt  p,  Re  '  ’ 


For  constant  viscosity,  the  divergence  of  the  dimensional  viscous  stress  tensor  can  be  written  as: 


V  •  r  =  -fiV  x  qj  +  AV  •  u 


(11) 


so  that 


V  x  (V  •  r)  =  -fN  x(Vxa) 
=fiV2d) 


(12) 


It  follows  that  Eq.  (10)  can  be  rewritten  as: 


doo  „  ,  .  Vo  Du  1 

+  Vx  oxu  = — — x  —  + - V2a) 

9t  p0  Dt  p0  Re 


(13) 


The  derivation  of  Eq.  (13)  completes  the  construction  of  the  physical  model;  a  summary  of  the 
governing  equations  in  given  in  Table  2.  For  clarity  and  convenience,  the  subscripts  o,  1  are 
omitted  from  Table  2;  tins  simplified  convention  will  also  be  adopted  in  the  remainder  of  this 
manuscript 


Table  2 

Summary  of  Model  Equations 


Helmholtz  Decomposition: 


u  =  V<f>  +  Vxys 
V2y  =  -co 


(2.1) 


i 


Continuity: 

< 

V20  =  -  — v2t 

y  p[  y  dt  Pe 

— —  =  -P\u*ndA  +  —  fVT  *ndA 
y  dt  J  Pei 

(2.2) 

State: 

P-pT 

(2.3) 

Energy: 

DT  _y~XdP  _  1  7,r 

Dt  y  dt  Pe 

(2.4) 

Vorticity  Transport: 

d(t)  „  ,  .  Vp  Du  1  r-,2 

-r-+  V  X  G)XU  = - — X - + - V  CO 

dt  p  Dt  p  Re 

(2.5) 

Normalization 

Despite  the  relative  simplicity  of  the  geometry  of  the  thermoacoustic  stack,  there  exists  a 
multitude  of  appropriate  conventions  which  can  be  adopted  in  the  nonnalization  of  the  governing 
equations.  Thus,  it  is  worthwhile  providing  clear  descriptions  of  normalizing  parameters,  and  of 
the  meaning  of  the  resulting  dimensionless  groups.  Below,  we  shall  choose  the  plate  separation 

distance,  H ,  as  reference  lengthscale  and  the  quantity  QH  as  reference  velocity  scale.  This 
convention  naturally  leads  to  the  following  definition  of  the  Reynolds  number: 


Re* 


(14) 


We  shall  refer  to  Re*  as  the  "Stokes  layer  Reynolds  number"  since  it  effectively  relates  the 
thickness  of  the  Stokes  layer, 


8  =*  6.4. 


(15) 


to  the  plate  spacing,  H ,  through: 


8_ 

H 


(16) 


Other  useful  conventions  for  the  Reynolds  number  include:  (a)  the  "acoustic  amplitude  Reynolds 
number". 


Rea» 
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(17) 


which  is  a  dynamic  Reynolds  number  based  on  the  acoustic  wave's  velocity  amplitude,  a,  the 
thickness  of  the  Stokes  layer,  and  the  kinematic  viscosity;  and  (b)  the  "particle  displacement 
Reynolds  number," 


Re, 


a 

M 


(18) 


which  characterizes  the  ratio  of  a  characteristic  particle  displacement  during  a  wave  cycle  to  the 
plate  separation  distance.  It  is  interesting  to  note  that  these  three  Reynolds  numbers  are  related 

through:  Rea2/  Rea  s  4 Re,2. 

Following  the  above  discussion,  the  physical  operating  conditions  of  a  thermoacoustic  stack 
are  determined  in  terms  of  the  following  set  of  dimensionless  parameters:  (a)  one  kinematic 
Reynolds  number,  i.e.  Rea  or  Rcp;  (b)  the  dynamic  Reynolds  number,  Rea  or,  equivalently,  the 
acoustic  drive  ratio,  R,  (c)  either  the  Prandtl  number,  Pr,  or  the  Pec  let  number,  Pe\  (d)  a  plate 
thickness  parameter,  e.g.  the  channel  width  to  plate  spacing  ratio,  h  /  H;  (e)  a  plate  length 
parameter,  e.g.  the  plate  length  to  spacing  ratio,  /  /  H ;  and  (0  the  location  of  the  thermoacoustic 
stack  with  respect  to  the  driving  acoustic  wave,  e.g.  in  terms  of  a  dimensionless  wavenumber  icc. 
Thus,  even  if  the  above  simplified  model  is  adopted,  a  high-dimensional  parameter  space  is  still 
needed  to  describe  the  relevant  dynamics. 


Computational  Methodology 

Simulation  of  the  governing  model  equations  is  performed  using  a  finite  difference 
methodology.  To  this  end,  the  computational  domain  is  divided  using  a  square  rectangular  grid 
with  mesh  size  Ax  =  Ay  (Fig.  3).  All  spatial  derivatives  are  discretized  using  the  standard  second- 
order  centered  differences,  while  time  derivatives  are  treated  using  the  third-order  Adams- 
Bashforth  integration  scheme. 


Periodicity  b.c.'s 


Plate 


Computational  Domain 


Grid 


Matching  to  acoustic  wave 


Figure  3.  Schematic  illustration  of  the  computational  domain. 
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Assuming  all  the  relevant  quantities  known  at  time  t  =  nAt,  the  solution  is  advanced  in  time  as 
summarized  below: 

(a)  Integrate  numerically  Eq.  (2.2b)  in  order  to  determine  the  thermodynamic  pressure  at  the  new 
time  step: 


pn+l  _  pn 

At 


v  J= o 


pn-, ju:-j.ndA+±j . n 


clA 


(19) 


(b)  Integrate  the  energy  equation  in  order  to  determine  the  new  temperature  distribution  within  the 
domain: 


r»n+I 


At 


-XA 


i= o 


U 


■  |  y-1  dP 

ypn~j  dt 


1 


+— V,T 
Pe 


2-T-/I-; 


(20) 


(c)  Using  the  equation  of  state,  Eq.  (2.3),  determine  the  new  value  of  the  density  distribution: 

n/*+ 1 

p"+!  (2D 

(d)  Based  on  the  new  values  of  pressure,  temperature  and  density,  compute  the  new  rate  of 
pressure  change  using  Eq.  (2.2b): 
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(e)  Solve  the  following  Neumann  problem  to  determine  the  new  potential  distribution  with  the 
computational  domain: 
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(f)  Integrate  the  vorticity  transport  equation  in  order  to  determine  the  new  vorticity  distribution  in 
die  interior  of  the  computational  domain: 
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(g)  Determine  the  new  streamfunction  distribution  within  the  computational  domain  and  the 
vorticity  distribution  along  solid  boundaries: 


[  y/**1  =0  on  all  boundaries 


(25) 


ACCOMPLISHMENTS 

Numerical  modelling  activities  conducted  during  the  first  phase  of  this  collaborative  research 
effort  resulted  in: 

(1)  development  of  the  simplified  physical  model  equations  summarized  above; 

(2)  formulation  of  a  computational  methodology  for  the  simulation  of  the  governing  equations; 
and, 

(3)  construction  of  computer  codes  reflecting  the  results  of  (1-2)  above,  and  vectorization  and 
optimization  of  one  version  of  the  code  for  the  CRAY-C90  at  the  Pittsburgh 
Supercomputer  Center. 

Implementation  of  the  numerical  schemes  have  so  far  focused  on  a  visualization  of  the 
flowfield  in  the  neighborhood  of  the  thermoacoustic  stack.  Samples  of  computations  performed 
are  presented  as  Figs.  4-6,  which  depict  instantaneous  streamfunction  distributions  for  different 
stack  locations  within  the  resonance  tube.  Detailed  analysis  of  the  results  of  these  computations 
will  be  given  elsewhere  [5], 

Future  plans  include: 

(1)  application  of  the  computational  schemes  to  analyze  the  impact  of  non-linear  dynamics  on 
flowfield  behavior ; 

(2)  analysis  of  computed  results  in  order  to  evaluate  one-dimensional  theories  and  to  derive,  if 
necessary,  collections  to  these  theories  in  order  to  account  for  neglected  dynamics; 

(3)  extension  of  the  physical  model  and  of  the  computational  schemes  in  order  to  (i)  account  for 
heating  loads,  and  (ii)  accommodate  a  more  detailed  coupling  between  flowfield  dynamics  in 
the  neighborhood  of  the  stack  and  the  driving  acoustic  wave;  and, 

(4)  validation  of  the  numerical  schemes  by  comparing  computed  predictions  with  experimental 
results;  and, 

(5)  implementation  of  the  validated  schemes  to  characterize  the  performance  of  thermoacoustic 
refrigerators  under  a  wide  range  of  operating  conditions. 
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NOMENCLATURE 


Roman 

a 

A 


c 


c 


p 


D 

D_ 

Dt 


~~  +  u  •  V 
dt 


h 

H 

k 

l 

L 

M 


n 


P 

P 


Pe  s  QH2  /  a 


Prs  via 
q  =  -kVT 
R 


Rea  =  2a/  t/vO 

R  cssQH2/v 

Rep=alQH 

t 


U  =  (u,v) 
V 


acoustic  wave  velocity  amplitude 

surface 

sound  speed 

specific  heat  at  constant  pressure 
computational  domain 

material  derivative 

channel  height 
plate  spacing 
thermal  conductivity 
plate  length 

length  of  computational  domain 
Mach  number 
outer  normal 
pressure 

thermodynamic  pressure 
Peclet  number 
Prandtl  number 
conduction  heat  flux 
drive  ratio 

acoustic  amplitude  Reynolds  number 

Stokes  layer  Reynolds  number 

particle  displacement  Reynolds  number 
time 

flow  velocity 
volume 
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thermal  diffusivity 


Greek 
a  =k  /  pc  p 

fi 

r 

5 

<p 

<p 

k  =  Cl  /  c 
X 

p 

v 

P 

i> 

CO 

¥ 

Ax,  Ay 
Cl 

Subscripts 

co 

o,l 

a 

h 

Superscripts 

n 


constant 

specific  heat  ratio 
thickness  of  the  Stokes  layer 
velocity  potential 
viscous  dissipation  function 
wavenumber 

second  coefficient  of  viscosity 
dynamic  viscosity 
kinematic  viscosity 
density 

shear  stress  tensor 
voracity 
streamf unction 
grid  size 

angular  frequency 


vortical  component 
zeroth  or  first  order 
acoustic 
discrete 


dimensional  quantity 
time  level 
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FIG  4:  Instantaneous  streamfunction  distribution  in  the  neighbourhood  of  a  thermoacoustic  stack  at 

(a)  0=  19.634,(b)  9  =  20.42,  (c)  9=  21.205,  (d)  9=  21.99,  (e)  9  =  22.77,  (0  0  =  23.56,  (g)  0=  24.34,  (fa)  0  = 

25.132.  The  stack  has  the  same  dimensions  as  stack  no  1.  as  used  by  Atchley  [6],  (L/H  =  4.727  h/H=0.727  )  and 

is  located  at  kx=rt/2  .  The  drive  ratio  R=  0.13%;R6a=  0.8879;Rep  =  0.0096;R6g  =  2132.  The  arrows  indicate 
the  magnitude  and  the  direction  of  the  acoustic  velocity  at  the  center  of  the  stack. 


HG  5:  Instantaneous  streamfunction  distribution  in  the  neighbourhood  of  a  thermoacoustic  stack  at 

(a)  0=  19.634,(b)  0=  20.42,  (c)  0  =  21.205,  (d)  0=  21.99,  (e)  0  =  22.77,  (0  0  =  23.56,  (g)  0  =  24.34,  (h)  0  = 

25.132.  The  stack  has  the  same  dimensions  as  stack  no  1.  as  used  by  Atchley  [  6  ],  (L/H  =  4.727  h/H=0.727  )  and 

is  located  at  kx=3tc/4  .  The  drive  ratio  R=  0.13%  ;  Rea=  0.8879  Rep  =  0.0096  ;  Re5  =  2132  .  The  arrows 
indicate  the  magnitude  and  the  direction  of  the  acoustic  velocity  at  the  center  of  the  stack 


X 


£1^6:  Instantaneous  streamfunction  distribution  in  the  neighbourhood  of  a  thermoacoustic  stack  at 

!£;634'<b>  0=  2042’  W  0=  2 1-205,  (d)  0=  21.99,  (e)  9  =  22.77,  (0  0  =  23.56,  (g)  0=  24.34,  (h)  0  = 
25. 132.  The  stack  has  the  same  dimensions  as  stack  no  1.  as  used  by  Atchley  [  6  ],  (L/H  =  4.727  h/H=0 .727  )  arid 

is  located  at  kx=n  .  The  drive  ratio  R=  0.13%  ;  Rea=  0.8879  Rep  =  0.0096  ;  Rt6  =  2132  .  The  arrows 
indicate  the  magnitude  and  the  direction  of  the  acoustic  velocity  at  the  matching  surfaces. 
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